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ABSTRACT 

We study the effects of anisotropic thermal conduction on low-cohisionahty, astrophys- 
ical plasmas using two and three-dimensional magnetohydrodynamic simulations. For 
weak magnetic fields, dilute plasmas are buoyantly unstable for either sign of the tem- 
perature gradient: the heat-flux-driven buoyancy instability (HBI) operates when the 
temperature increases with radius while the magnetothermal instability (MTI) operates 
in the opposite limit. In contrast to previous results, we show that, in the presence of 
a sustained temperature gradient, the MTI drives strong turbulence and operates as 
an efficient magnetic dynamo (akin to standard, adiabatic convection). Together, the 
turbulent and magnetic energies contribute up to '^10% of the pressure support in the 
plasma. In addition, the MTI drives a large convective heat flux, ~ 1.5% xpc^. These 
findings are robust even in the presence of an external source of strong turbulence. Our 
results on the nonlinear saturation of the HBI are consistent with previous studies but 
we explain physically why the HBI saturates quiescently by re-orienting the magnetic 
field (suppressing the conductive heat flux through the plasma), while the MTI satu- 
rates by generating sustained turbulence. We also systematically study how an external 
source of turbulence affects the saturation of the HBI: such turbulence can disrupt the 
HBI only on scales where the shearing rate of the turbulence is faster than the growth 
rate of the HBI. In particular, our results provide a simple mapping between the level of 
turbulence in a plasma and the effective isotropic thermal conductivity. We discuss the 
astrophysical implications of these findings, with a particular focus on the intracluster 
medium of galaxy clusters. 

Key words: conduction, convection, instabilities, MHD, turbulence — galaxies: clus- 
ters: general 



1 INTRODUCTION 

The thermodynamics of a plasma can have dramatic and 
sometimes unexpected implications for its dynamical evolu- 
tion. For example, thermal conduction can reduce the accre- 
tion rate in spherical accretion flows by as much as two to 



three orders of magnitude relative to the Bondi value ( [John- 
son k Quataert 2007; Shcherbakov Baganoff 201Qt. More 



relevant to this paper, Balbus (2000) and Quataert 



(2008) 



demonstrated that the convective dynamics of conducting 



* E-mail:mkmcc@astro. berkeley.edu 



plasmas are completely different from those of an adiabatic 
fluid. This paper focuses on the nonlinear evolution and sat- 
uration of this convection. 

When anisotropic conduction is rapid compared to the 
dynamical response of a plasma, the temperature gradient, 
rather than the entropy gradient, determines the plasma's 
convective stability. The convective instability in this limit 
is known as the heat-flux- driven buoyancy instability (HBI) 
when the temperature increases with height {g • VT < 0) or 
the magnetothermal instability (MTI) when the temperature 
decreases with height {g • VT > 0). We summarize the lin- 
ear physics of these instabilities in section [2l [Parrish Stone] 
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|2005';2007') andlParrish & Quataertl f 2008) studied the non- 



linear development of the MTI and HBI using numerical sim- 
ulations. These instabilities couple the magnetic structure of 
the plasma to its thermal properties and potentially have im- 
portant implications for galaxy cluste rs (Parrish et al.||2008 



2009; BogdanovTc et al.||2009; Sharma et al. |2009a| |Parris 



et al..2010, Bogdanovic et al.,,2010, ^Ruszkowski Oh 2010 



hot accretion flows onto compact objects ( Sharma et al 



2008 



and the interiors and surface layers of white dwarfs and neu- 
tron stars ( ,Chang et a l. 2 010^) . 

In this paper, we revisit the nonlinear behavior of the 
HBI and MTI. We focus on the physics of their saturation, 
but also include a lengthy discussion of possible astrophys- 
ical implications in section |6] Our analysis is idealized (we 
use a plane-parallel approximation and neglect radiative cool- 
ing), but we are able to understand the nonlinear behavior of 
the HBI and MTI, and therefore their astrophysical implica- 
tions, more thoroughly than in previous papers. Our results 
for the saturation of the HBI are similar to those of Parr ish fcl 



Quataert ( ,2008) , but with an improved understanding of the 
saturation mechanism. Our results for the MTI, however, dif- 
fer from previous results, significantly changing the predicted 
astrophysical implications of the MTI. We show in section [4?2] 
that this is because the development of the MTI in many pre- 
vious simulations has been hindered by the finite size of the 
simulation domain. 

The structure of this paper is as follows. We describe the 
linear physics of the MTI and HBI in section [2] We describe 
our computational setup in section |3] and the results of our 
numerical simulations in section [4] In order to better under- 
stand how the saturation of the HBI and MTI may change 
in a more realistic astrophysical environment, we study the 
interaction between these instabilities and an external source 
of turbulence or fluid motion (section [5]). Finally, in section [g] 
we summarize our results and discuss the astrophysical im- 
plications of our work. 



2 BACKGROUND 

2.1 Equations and Assumptions 

We assume that the plasma is an ideal gas with an adiabatic 
index 7 = 5/3 and model it using the magneto- hydrodynamic 
(MHD) equations, neglecting all dissipative processes except 
thermal conduction. In Cartesian coordinates, the equations 
for the conservation of mass, momentum and magnetic flux, 
and for the evolution of internal energy are 

dp 



dt 



+ V - (pi 



0, 



d_ 
dt 



pT 



OB 

ds 



dt 



pv ^ V -\- 
V X {v X B 

^ Qcond' 



(la) 
(lb) 

(Ic) 
(Id) 



where p is the mass density, v is the fluid velocity, denotes a 
tensor product, P is the pressure, B is the magnetic field, I is 
the unit matrix, g is the gravitational field, / is an externally 
imposed force, T is the temperature. 



kB 



7 — 1 rriH 



In 



(2) 



is the entropy per unit mass, d/dt = d/dt + v - V is the 
Lagrangian time derivative, and Qcond conductive heat 

flux. 

We ignore the ion contribution to the conductive heat 
flux, which is smaller than the electron contribution by a fac- 
tor of (mi/me)^/^ ^ 42. We assume that the electrons have 
mean free paths much longer than their gyro-radii (as is the 
case in the intracluster medium (ICM) of galaxy clusters); the 
electrons therefore move almost entirely along magnetic field 
lines. Consequently, the thermal conductivity of the plasma 
is strongly anisotropic. The conductive heat flux in this limit 
becomes 



Qco 



-Keb(b- VT), 



(3) 



where h — B / B \s ?l unit vector in the direction of the mag- 
netic field and is the thermal conductivity of free electrons 
( Braginskii||1965|. W hile depends sensitively on temper- 
ature (Spitzer 1962), we take it to be constant in our cal- 
culations to simplify the interpretation of our results. This 
approximation does not affect our conclusions, because the 
physics of the buoyancy instabilities is independent of the con- 
ductivity in the limit that the thermal diffusion time across 
the spatial scales of interest is short compared to the dynam- 
ical time. 



2.2 The Physics of Buoyancy Instabilities in Dilute 
Plasmas 



The dissipative term in equation (Id) shows that fluid dis- 



placements in a plasma are not in general adiabatic. The stan- 
dard analysis of buoyancy instabilities therefore does not ap- 
ply, and the convective and mixing properties of a conducting 
plasma can be very different from those of an adiabatic fluid. 
In this paper, we focus on the limit in which thermal conduc- 
tion is much faster than the dynamical time in the plasma. 
This applies on scales < 7(Aif)^/^, where A is the electron 
mean free path and H is the plasma scale height; this "rapid 
conduction limit" encompasses most scales of interest for the 
ICM. In this limit, the magnitude of the temperature gra- 
dient and the local orientation of the magnetic field control 
the convective stability of the plasma ( |Balbus|2000 Quataert 
20081. 



Although a single dispersion relation describes the lin- 
ear stability of plasmas in this limit, it is easiest to un- 
derstand the physics by separately considering cases where 
the temperature increases or decreases with height. |Balbus| 
(2000) first considered the case where the temperature de- 



creases with height and identified the magnetothermal insta- 
bility, or MTI. We show a schematic of this instability in the 
top row of Figure [l] The first panel of this figure shows a 
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Figure 1. Illustration of the linear development of buoyancy in- 
stabilities in dilute plasmas from two-dimensional numerical simu- 
lations. Color shows the temperature, increasing from blue to red, 
while black lines trace the magnetic field. Top Row: Quasi-linear 
evolution of the MTI. Left panel: Initial equilibrium state, with 
bz = and the conductive heat flux Q = everywhere. Right 
panel: The plasma at t = 5 tbuoy given an initial perturbation with 
k = 27t/L X. Flux freezing and rapid conduction along fleld lines 
ensure that Lagrangian fluid displacements are nearly isothermal. 
If the initial state has a positive temperature gradient, upwardly 
displaced fluid elements are warmer than their surroundings; they 
will expand and continue to rise. Similarly, downwardly displaced 
fluid elements are cooler than their surroundings and will sink. 
Bottom Row: Quasi-linear evolution of the HBI. The initial equi- 
librium state has bz = 1 and a net downward conductive heat flux, 
with V • Q = everywhere. Left panel: The plasma at t = 5 tbuoy 
given an initial perturbation with kx = kz- The perturbation alters 
the geometry of the magnetic fleld, so that fleld lines are not paral- 
lel and there can be conductive heating and cooling of the plasma. 
Right panel: Temperature difference at t = 5 tbuoy relative to the 
initial condition. Magnetic fleld lines converge in upwardly dis- 
placed fluid elements, leading to an increased temperature. These 
fluid elements will then expand and continue to rise. Similarly, 
downwardly displaced fluid elements are conductively cooled, con- 
tract and continue to sink. 



plasma in hydrostatic and thermal equilibrium, with a weak 
horizontal magnetic field. We apply a small, plane wave per- 
turbation and, as the plasma evolves, the field lines follow 
the fluid displacements. Efficient conduction along these field 
lines keeps the displacements isothermal. Since the tempera- 
ture falls with height, upwardly displaced fluid elements are 
warmer than their new surroundings; they expand and con- 
tinue to rise. Similarly, downwardly displaced fluid elements 
sink. An order of magnitude calculation shows that displace- 



ments grow exponentially, with the growth rate (or e-folding 
rate) p ^ \g d\nT/ dz 



-1/2 



Quataert ( 2008 ) investigated the limit in which the tem- 



perature of the plasma increases with height. The instability 
in this case is known as the heat-flux-driven buoyancy insta- 
bility, or HBI. We sketch the growth of the HBI in the bottom 
row of Figure [l] Here, the initial equilibrium state has vertical 
magnetic field lines and a constant heat flux; the latter is re- 
quired for the plasma to be in thermal equilibrium. However, 
perturbations to the magnetic field divert the heat flux and 
conductively heat or cool pockets of the plasma. If the temper- 
ature increases with height, upwardly displaced fluid elements 
become warmer than their surroundings and therefore expe- 
rience a destabilizing buoyant response. These perturbations 
grow with the same growth rate \g d\nT/dz\~'^^'^ . 

iQuataerT ( |2008| performed a WKB analysis on equa- 
tions (la)-(ld) in the Boussinesq limit and obtained the fol- 



lowing dispersion relation for plasma in a constant gravita- 
tional field g = —gz, threaded by a constant magnetic field 
with any orientation in the x — z plane: 



dz 



[(2b^-l)(l. 



ki 



'^bxbzkxkz 



(4) 



Here, p = —iu is the local growth rate of the mode (the time 
dependence of the perturbation is e^*), uja = k ■ va is the 
Alfven crossing frequency. 



7 — 1 /iTTlH ds 

7 kB ^ dz 



(5) 



describes the buoyant response of an adiabatic plasma, k 
k/k is the direction of the wave vector, and 

7-1 KeT ff^ \2 

is inversely proportional to the conduction time across the 
wavelength of the mode. It is convenient to identify cjbuoy = 
\g dhiT / dz\^^'^ as a characteristic frequency for the buoyancy 



(6) 



buoy 



and ta 



N- 



instabilities and we will also use tbuoy 
in our analysis. 

In the limit that conduction is rapid compared to any 
dynamical response (cj^, ^ cja, cjbuoy), equation Q sim- 
plifies tc[^ 

(/ + uj\) = sgn{dT/dz) cjbuoy x 

(7) 



[(2^?-l)(l 



kt 



^bxbzkxk 



As we stated in section [23] the growth rate in the rapid con- 
duction limit is independent of the thermal conductivity. 

Equation ^ shows that magnetic tension can suppress 
the MTI and HBI; if cj^ > <^buoy5 must be negative, and 
the plasma is stable to small perturbations. In this paper, 
we focus on the relatively weak field limit in which magnetic 

^ In going from the cubic equation Q to the quadratic ([7|, we have 
ignored a solution for p. This solution is exponentially damped on 
the conduction timescale and is not relevant to our analysis. 
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tension does not suppress the instabilities, and we take uja ^ 
cjbuoy in equation ([7|. The dominant role of the magnetic 
field in our analysis is to enforce anisotropic electron heat 
transport. 

For any magnetic field direction 6, the term in square 
brackets in equation ([7| can be positive or negative; thus, 
there are always linearly unstable modes, irrespective of the 
thermal state of the plasma (excluding the singular case of 

dT/dz^^). 

We finally note that, although we neglected cooling in 
term to equation (|7| which is in- 



equation (Id), it adds 



significant when the cooling time is much longer than uj^ ^ 
( Balbus Reynolds|[2008l|201Q| . This is not necessarily the 
case near the centers of cool-core clusters, and we plan to 
explore the combined effects of cooling and buoyancy insta- 
bilities in future work. For now, we assume that the buoyancy 
instabilities develop independently of cooling; this allows us 
to understand the nonlinear development and saturation of 
the instabilities themselves and therefore to assess their pos- 
sible implications for the thermal balance of the plasma. 



3 NUMERICAL METHOD 

3.1 Problem Setup and Integration 

We consider the evolution of a volume of plasma initially in 
hydrostatic and thermal equilibrium, but subject to either the 
HBI or MTI. We seed our simulations with Gaussian-random 
velocity perturbations with a flat spatial power spectrum and 
a standard deviation of 10~^ Cg. The small amplitude of these 
initial perturbations ensures that the instabilities start out in 
a linear phase and permits us to compare our results with the 
predictions of equation ([7|. Astrophysical perturbations are 
unlikely to be this subsonic, however, and the instabilities in 
our simulations take much longer to saturate then one would 
expect from the larger perturbations found in a more realistic 
scenario. We return to this point in section |6] 

We highhght the distinctness of the HBI and MTI 
from adiabatic convection by choosing initial conditions with 
ds/dz > whenever possible, so that the plasma would be 
absolutely stable if it were adiabatic. Our results do not rely 
critically on this choice, however, because the plasma is not 
adiabatic and its evolution is independent of its entropy gra- 
dient to lowest order in c^buoy/^^- Our results are much easier 
to interpret when cjbuoy doesn't vary across the simulation do- 
main, and we prioritize this constraint on the initial condition 
over the sign of its entropy gradient (as we discuss in more 
detail below). 



more computationally expensive than adiabatic MHD calcu- 
lations. We draw most of our conclusions from simulations 
performed on uniform Cartesian grids of (64)^ and (128)^. 
We also performed a large number of two-dimensional simu- 
lations on grids of (64)^ (128)^ and (256)^ We found that 
the kinetic energy generated in our local HBI simulations con- 
verges by a resolution of (64)^; our global HBI simulations 
require roughly 60 grid cells per scale height to give a con- 
verged kinetic energy. The MTI is somewhat more sensitive 
to resolution; we tested the convergence of these simulations 
by comparing the results of an identical simulation performed 
at (128)^ and (256)^ 

We perform our calculations in the plane-parallel approx- 
imation with uniform gravity g = —gz. We fix the temper- 
ature at the upper and lower boundaries of our computa- 
tional domain, and we extrapolate the pressure into the upper 
and lower ghost cells to ensure that hydrostatic equilibrium 
holds at the boundary. For all other plasma variables, we ap- 
ply reflecting boundary conditions in the direction parallel 
to gravity and periodic boundary conditions in the orthogo- 
nal directions. As we discuss in section |4] the choice of fixed 
temperatures at the upper and lower boundaries has an im- 
portant effect on the non-linear evolution. Simulations with 
Neumann boundary conditions in which the temperature at 
the boundaries is free to adjust would give somewhat dif- 
ferent results (see, e.g., Parrish Stone|2005 ). Our choice of 
Dirichlet boundary conditions is largely motivated by the fact 
that many galaxy clusters in the local universe are observed 



to have non-negligible temperature gradients ( [Piffaretti et al. 
'20051 



We perform both local (with the size of the simulation 
domain L much smaller than the scale height H) and global 
(L > H) simulations of the MTI and HBI. The local simu- 
lations separate the development of the instability from any 
large-scale response of the plasma, allowing us to study the 
dynamics in great detail. However, because the response of 
the plasma on larger scales can influence the nonlinear evo- 
lution and saturation of the instabilities, we also carry out 
global simulations. 



3.2 Local Simulations 



In sections |3j|5] we work in units with = firrip = 1. As 
noted previously, we restrict our analysis to plasmas with 
rapid conduction; we find that setting /€e = 10 x pcjbuoy^^ 
puts our local simulations safely in this limit. This corre- 
sponds to a thermal diffusion time across the box of ^ 
We solve equations 0-}!^, along with a conservative 0-1 ^wy We initiahze our local HBI simulations with a hnear 



form of equation (Id), using the conservative MHD code 
Athena ( [Gardiner Stone] [2Q08| [Stone et a'r]|2008t with 
the anisotropic conduction algorithm described in [Parrish[ 
k Stone[ (|2005) and Sharma k Hammett (2007). In partic- 



buoy ■ 

temperature gradient: 



ular, we use the monotonized central difference limiter on 
transverse heat fluxes to ensure stability. This conduction 
algorithm is sub-cycled with respect to the main integrator 
with a time step At oc (Ax)^; our simulations are therefore 



T{z) = To (1 + z/H) , 
piz) = po{l + z/H)-^, 
P(z) = p{z) T{z) . 



(8b) 



We choose to set po = 1 and g = 1 = 2To/H. Unless other- 
wise noted, we take H — 2 (To = 1) in our local simulations 



© 0000 RAS, MNRAS 000, 000-000 



Nonlinear Saturation of Buoyancy Instabilities 5 



and we evolve a volume of plasma from z = to z 
i.e., over ^ 5% of a scale-height. 



0.1, 



atmosphere defined by equation (10b) does not have a sin- 



We use the setup of Parrish & Stone ( 2005 ) for our local 
MTI simulations: 



Tiz) = To (1 - z/H) , 
p{z) = po (1 - z/Hf , 
P{z) = p{z) T{z), 



(9a) 
(9b) 
(9c) 



with g = 1 and H — ?>. The MTI induces large vertical dis- 
placements in the plasma, and our reflecting boundary condi- 
tions clearly influence its evolution; we attempt to minimize 
the effects of the boundaries by sandwiching the unstable vol- 
ume of p lasma between two buoyantly neutral layers. parrish| 
&: Stone| (f2005 ) describe this setup in more detail]^ 

Both of the above atmospheres have positive entropy 
gradients and therefore would be stable in the absence of 
anisotropic conduction. The results we describe in this pa- 
per are entirely due to non-adiabatic processes. 

We show in section |4.1| that local and global HBI sim- 
ulations give very similar results; we therefore use the sim- 
pler, local simulations for most of our analysis of the HBI. 
By contrast, the large-scale vertical motions induced by the 
MTI make its evolution inherently global; as we describe in 
section [X2] local simulations do not give a converged result 
independent of L/ H. We therefore study the saturated state 
of the MTI using global simulations, which require an alter- 
nate set of initial conditions. 



3.3 Global Simulations 

The physical properties of the plasma in a global simulation 
can vary by an order of magnitude or more across the simula- 
tion domain. This complicates our study of the HBI and MTI 
because the instabilities can be in different stages of evolu- 
tion at different spatial locations. Similarly, the kinetic and 
magnetic energies generated by the instabilities in a global 
simulation can be functions of height, obscuring their depen- 
dence on other parameters. While one must confront these 
problems when studying buoyancy instabilities in an astro- 
physical context, we avoid them by choosing initial conditions 
with a buoyancy time that is constant with height. We specify 



g{z) go e 
T{z) exp 



-z/S 



90 



(10a) 
(10b) 



The positive and negative signs in equation (10b) produce 



atmospheres unstable to the HBI and MTI, respectively. Note 
that T(z = 0) = 1 as in our local simulations. We take = 1, 
S = 3 and c^^uoy — 1/2 and numerically solve for p{z) so 
that the initial atmosphere is in hydrostatic equilibrium. The 



2 Note that, although [Parrish Stone] ( |2005| > describe these extra 
layers as buoyantly stable, the effect of the isotropic conductivity 
is to make them buoyantly neutral. 



gle, well-defined scale height. We therefore define H so that 
L/H — In (Tmax/^lnin); thls definition for H is analogous to 
that used in our local simulations and L/H reflects the total 
free energy available to the instabilities. Because if is a de- 
fined, rather than fundamental, property of the atmosphere, 
it is not independent of the size of the simulation domain. 
We carry out simulations with domain sizes L = 0.5 and 2.0, 
which have L/H — 0.27 and 1.4, respectively. 

Simulations with a larger simulation domain size inher- 
ently permit larger vertical displacements, and the boundary 
conditions do not influence the evolution of the plasma as 
strongly as they do in a local simulation. We find that the 
neutrally stable layers described in section [3^ don't alter the 
results of our global simulations, and thus we do not include 
them in our setup. 

We chose the conductivity in our local simulations so 
that = lOpcJbuoyi^^, but the corresponding constraint on 
the time step becomes impractical for our global simulations. 
Instead, we adjust so that the ratio Ke/L is the same as it is 
in the local simulations. This keeps the ratio of the conduction 
time to the sound crossing time constant, and is appropriate 
if the turbulence driven by the MTI or HBI reaches a terminal 
speed less than or of order the sound speed. We carried out 
simulations with different values of the conductivity and ver- 
ified that our results are insensitive to factor of few changes 

in K,e- 

The atmosphere we use for our global MTI simulations 
has a negative entropy gradient. One might worry that the 
results of these simulations — which aim to focus on the MTI — 
would be biased by the presence of adiabatic convection. This 
is, however, only a pedagogical inconvenience; cj,^ > on 
all relevant scales in the simulation, so the adiabatic limit to 
equation (|4| is not important, and the evolution of the plasma 
is nearly independent of its entropy gradient. To confirm this 
we carried out L — 0.27 H simulations with both sets of initial 
conditions; they give very similar results. Additionally, we 
performed one simulation of adiabatic convection using this 
setup, and the behavior of this simulation is entirely different 
from the MTI, especially at late times. 



3.4 Turbulence 

The pure HBI/MTI simulations described in the previous 
sections are physically very instructive but astrophysically 
somewhat idealized. In order to better understand the as- 
trophysical role of the HBI and MTI, we also study their 
interaction with other sources of turbulence and fluid mo- 
tion. We assume that these take the form of isotropic tur- 
bulence; we include such turbulence via the externally im- 
posed force field / in equation (lb). Following Lemaster 
Stone| ( |2008| , we compute / in momentum space from scales 
ko = {4,6,8} X 27r/L, down to kma^ = 2/co, with an injected 



energy spectrum E^^^'^ oc k~'^. We then randomize the phases, 
perform a Helmholtz decomposition of the field, discard the 
compressive component, transform the field into configura- 
tion space, and normalize it to a specified energy injection 
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Table 1. Parameters for the HBI simulations (§ |4.1 



Name 


D 


res 


L/H 




hi 


2 


64 


0.05 


7.07 


h2 


3 


64 


0.05 


7.07 


h3 


3 


128 


0.75 


0.47 


h4 


3 


128 


1.40 


0.35 



All simulations are performed on square Cartesian grids of size 
L. D is the dimensionality of the simulation, res is the number 
of grid cells along a side, and k, is the conductivity (in units of 
k^/lim X pc^buoy^^)- All of these simulations are local (eq. [sj, 
except for the one with L/H = 1.4, which uses the global setup 
(eq. |10| >. We initialized all of these simulations with weak horizontal 
magnetic fields (B/ViTT = 10-^). 



rate. We find that the turbulence sets up a nonlinear cascade 
that is not very sensitive to either the driving scale or the 
injected spectrum of the turbulence. 

This prescription for turbulence is statistically uniform 
in both space and time and therefore provides a controlled 
environment in which to study the interaction between tur- 
bulence and buoyancy instabilities. This may not, however, 
be a good approximation to real astrophysical turbulence; we 
discuss the implications of our choice in sections |5] and |6] 



4 NONLINEAR SATURATION 
4.1 Saturation of the HBI 



Parrish & Quataert ( 2008 ) described the nonlinear saturation 
of the HBI, but they did not explicitly test the dependence 
of their results on the size of the computational domain. This 



dependence turns out to be crucial for the MTI (see § 4.2 ), but 
we show here that the size of the domain has little effect on the 
saturation of the HBI. Nonetheless, we describe the nonlinear 
behavior of the HBI in reasonable detail, expanding on the 
physical interpretation given in previous papers (although our 
HBI results do not differ significantly from those of previous 
authors). The saturation of the HBI is the simplest process 
we consider in this paper and serves as a useful comparison 
for our new results. 

We study the saturation of the HBI using 2D and 3D 
simulations spanning a range of domain sizes L/H. Table [l] 
lists all of the simulations presented in this section. 

Figure [2] shows snapshots of the evolution of temperature 
and magnetic field lines in a local, 2D HBI simulation. We 
chose this simulation to simplify the field-line visualization, 
but the results in Figure [2] apply equally to our local and 
global 3D simulations. We initialized this simulation in an 
unstable equilibrium state with ver tica l magnetic field lines 

we seed this initial 



3.1 



{hz — 1). As described in section 
condition with small velocity perturbations; the HBI causes 
these perturbations to grow in the first three panels of Fig- 
ure [2] The evolution becomes nonlinear in the third panel, 
when the velocity perturbations reach ^ 4% of the sound 



speed. Afterwards, the instability begins to saturate and the 
plasma slowly settles into a new equilibrium state. The last 
panel in Figure |2] shows that this saturated state is highly 
anisotropic: the magnetic field lines are almost entirely or- 
thogonal to gravity. Flux conservation implies that the fluid 
motions must also be anisotropic, with most of the kinetic en- 
ergy in horizontal motions at late times (see Fig. [s] discussed 
below). These horizontal motions are very subsonic: in all of 
our simulations, the velocities generated by the HBI are sig- 
nificantly less than 1% of the sound speed in the saturated 
state. 

Because the fluid velocities remain small, the linear dis- 
persion relation (eq. [?]) captures much of the evolution of the 
HBI, even at late times. For any magnetic fleld orientation, 
the fastest growing modes are the ones with k along the axis 
h X (b X g); these modes have the growth rate 

Pmax = IcJbuoy bz\, (11) 

which decreases as the field lines become horizontal. Addi- 
tionally, when bl < 1/2, only modes with /c^ > 1 — 4(6^ — bi) 
are unstable. Since the HBI saturates by making the field 
lines horizontal {bz 0), both the maximum growth rate 
of the instability and the volume of phase space for unsta- 
ble modes decrease as the HBI develops. This strongly limits 
the growth of the perturbations, and helps explain why the 
instability saturates relatively quiescently. 



As argued by Parrish & Quataert ( 2008 ) , the HBI satu- 
rates when its maximum growth rate pmax vanishes, so that 
there are no longer unstable modes. While this is clearly a 
sufficient condition for the plasma to reach a new stable equi- 
librium, it is by no means necessary: the instability could, e.g., 
saturate via nonlinear effects, but in practice this is not the 
case (at least for simulations without an additional source 
of turbulence). Equation (11) for pmax shows that the HBI 
could saturate by making either dT/dz or bz vanish; intu- 
itively, the HBI is powered by a conductive heat flux, which 
it must extinguish in order to stop growing. Erasing the tem- 
perature gradient might seem like the more natural satura- 
tion channel, since the conduction time across the domain is 
much shorter than the time it takes the HBI to develop and 
saturate. In an astrophysical setting, however, the large-scale 
temperature field is often controlled by cooling, accretion or 
other processes apart from the HBI. We therefore impose the 
overall temperature gradient in our simulations by fixing the 
temperature at the top and bottom of the domain, so that 
<^buoy is roughly independent of time and saturation requires 
bz=0. 

Since the HBI saturates by making the magnetic field 
lines horizontal, we take the bz ^ limit in equation ([7| to 
understand the late-time behavior of the plasma: 

UO = ±UJhuoy {1 - kzj . (12) 

We have assumed here that the magnetic field is weak enough 
for magnetic tension to be negligible on the scales of inter- 
est [^Equation (12) shows that the saturated state of the HBI 



Note that equation |l2l only strictly applies when the magnetic 
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Figure 2. Evolution of the HBI with an initially vertical magnetic field in a local, 2D simulation (simulation hi in Table [T|. Color shows 
temperature and black lines show magnetic field lines. A small velocity perturbation to the initial state seeds exponentially growing modes 
which dramatically reorient the magnetic field to be predominantly horizontal. The induced velocities are always highly subsonic and, after 
t ~ 20tbuoy5 cire also almost entirely horizontal. Once the plasma reaches its saturated state, it is buoyantly stable to vertical displacements. 
The plasma does not resist horizontal displacements, but the saturated state is nearly symmetric to these displacements and they do not 
change its character. 



is buoyantly stable, but that there is a family modes with 
kz — 1 which feel no restoring force. This simply reflects the 
fact that the plasma is stably stratified and resists vertical 
displacements. Displacements orthogonal to gravity are un- 
affected by buoyancy, however, and appear as zero-frequency 
modes in the dispersion relation. 

Figure |3] demonstrates this asymmetry between vertical 
and horizontal displacements in the late time evolution of the 
plasma. This figure shows the kinetic energy in horizontal and 
vertical motions as a function of time. During the initial, lin- 
ear growth of the instability (t < 10 tbuoy), buoyantly unstable 
fluid elements accelerate toward the stable equilibrium, and 
the kinetic energy is approximately evenly split among verti- 
cal and horizontal motions. As the instability saturates, how- 
ever, the plasma becomes buoyantly stable and traps the ver- 
tical motions in decaying oscillations (internal gravity waves). 
The horizontal motions keep going, however, and retain their 
kinetic energy for the duration of the simulation. This differ- 
ence in the response of the plasma to vertical and horizontal 
motions accounts for the anisotropy of the velocity field in 
the saturated state of the HBL 

Figure |4] shows the evolution of the rms magnetic field 
angle (left panel) and magnetic energy (right panel) in 3D 
HBI simulations for three different values of the size of the 
computational domain L relative to the scale height H. The 
zero frequency modes discussed above also dominate the evo- 
lution of the magnetic field at late times, after the motions be- 
come nonlinear (t > 10 tbuoy). The horizontal displacements 
stretch out the field lines, amplifying and reorienting them. 
Quantitatively, we expect that bz ^ \/ ^ (x. t~^, where A is 
a characteristic scale for the modes in the saturated state, 
^ is the magnitude of the horizontal displacements, and we 
have assumed that the velocity is constant with time. The 
left panel of Figure [4] shows that the dependence in the sim- 



ulations is quite close to this, with bz oc t~°"^^[^ Stretching 
the field lines in this manner amplifies the field strength by 
an amount 5B (x. ^] \i the velocity is constant with time, we 
expect oc t^. The right panel of Figure [Z] shows that this 
time dependence is approximately true for our two larger HBI 
simulations; the amplification is slightly slower in the very lo- 
cal calculation with L/H — 0.05. There is no indication that 
the magnetic field amplification has saturated at late times 
in the HBI simulations. We suspect that the amplification 
would continue until the magnetic and kinetic energy densi- 
ties reach approximate equipartition, but we would have to 
run the simulation for a very long time to verify this. 

One of the important results in Figure [4] is that the sat- 
urated state of the HBI is nearly independent of the size of 
the computational domain L/H. Since the late time evolution 
of the plasma is driven only by horizontal displacements, the 
key dynamics all occur at approximately the same height in 
the atmosphere. The saturation of the HBI is thus essentially 
local in nature and should not be sensitive to the global ther- 
mal state of the plasma or the details of the computational 
setup. 

The dramatic reorienting of the magnetic field caused by 
the HBI severely suppresses the conductive heat flux through 
the plasma. The conductive flux is proportional to (6^), which 
decreases in time oc (t/tbuoy)""*^"^- The saturated state of the 
HBI is also buoyantly stable and resists any vertical mixing 
of the plasma. As a result, the convective energy fluxes in our 
HBI simulations are very small, - 10~^ pel. The effect of the 
HBI therefore is to strongly insulate the plasma against both 
conductive and convective energy transport. This can dra- 
matically affect the thermal evolution of the plasma ( Parrish| 
et al.||2009||Bogdanovic et al.||2009| ). 



The fact that the growth rate of the HBI depends on the 
local orientation of the magnetic field, as well as the thermal 
structure of the plasma, makes it very different from adiabatic 



field is exactly horizontal. More generally, there will still be unstable 
modes with growth rate given by equation |11| this growth rate is 
very slow in the saturated state of the HBI, however, and these 
modes don't change the dynamics of the plasma. 



^ The slight difference relative to the simple predictions of fiux 
freezing given the velocity field in Figure [3] may be due to the finite 
resolution of our simulations, which prevents us from resolving the 
field line direction when 6^ < 10 x dz/L. 
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Figure 3. Evolution of the vertical and horizontal kinetic energy 
in a local, 2D HBI simulation (simulation hi in TableQ. The units 
are such that the thermal pressure P ^ 1 and the initial magnetic 
energy is S^/Stt = 10"-*^^. After a period of exponential growth 
in which the x and z motions are in approximate equipartition, 
the HBI saturates and the kinetic energy ceases to grow. At this 
point, the energy in the vertical motion is in the form of stable 
oscillations, which decay non-linearly. The horizontal motions are 
unhindered, however, and persist for the entire duration of the 
simulation. These horizontal motions are responsible for the asym- 
metry of the magnetic field shown in Figure^ 



convection. This dependence on the magnetic field structure 
provides a saturation channel in which the kinetic energies are 
very small compared to the thermal energy (e.g., in Fig. [s] 
pv^ /nkT ^ 10~^). Critically, these highly subsonic motions 
occur in simulations in which the boundary conditions allow 
for the presence of a sustained, order unity, temperature gra- 
dient {d\nT / dhi z ^ 1). In an adiabatic simulation, the anal- 
ogous sustained entropy gradient would generate convective 
motions with pv^ ^ nkT. This does not occur in an HBI- 
unstable plasma. Thus, although a plasma with a positive 
temperature gradient is in general buoyantly unstable, the ef- 
fect of the HBI is to peacefully stabilize the plasma within a 
few buoyancy times by suppressing the conductive heat flux 
through the plasma. The resulting, stably-stratified plasma 
then resists vertical mixing and, in the absence of strong ex- 
ternal forcing, we expect the fluid velocities and magnetic field 
lines to be primarily horizontal. In section |5] we perturb this 
state with externally driven, isotropic turbulence and test the 
strength of the stabilizing force. 



4.2 Saturation of the MTI 

Figure [5] shows the evolution of one of our local, 2D MTI 
simulations. As in the HBI simulation shown in Figure [2] we 



TTTT| 1 I I lllllj 1 I I MM 

- L/H = 1.4 - 

L/H = 0.75 - 

L/H = 0.05 




LUlL 



^ pTT| 1 I I lllllj 1 I I IIIU 



00 



cq 




LUlL 



10 

^/ ^buoy 



102 



10 

^/ ^buoy 



102 



Figure 4. Evolution of the orientation (left) and energy (right) of 
the magnetic field in 3D HBI simulations for three different values 
of the size of simulation domain relative to the temperature scale- 
height (L/H) (simulations h2-h4 in Table Units are such that 
the thermal pressure P is ~ 1. The results are nearly independent 
of size of the simulation domain. As described in § |4.1| the HBI 
saturates by shutting itself off; the linear exponential growth ends 
at f ~ 10 tbuoy 5 and most of the evolution of the magnetic field hap- 
pens afterward. This evolution is driven by the horizontal motions 
shown in Figure [3] which both amplify and reorient the magnetic 
field. After a brief period of exponential growth, the field amplifi- 
cation is roughly linear in time. By contrast, the field amplification 
by the MTI is exponential in time (see Fig, [s}. 



initialized this simulation in an unstable equilibrium state (a 
weak horizontal magnetic field) and seeded it with the small 
velocity perturbations described in section [3T] The MTI and 
HBI stem from very similar physics, and as a result have very 
similar linear dynamics. The nonlinear behavior of the two 
instabilities is entirely different, however. While the HBI sat- 
urates relatively quiescently by driving the plasma to a buoy- 
antly stable and highly anisotropic state, the MTI generates 
vigorous, sustained convection that tends to isotropize both 
the magnetic and velocity fields. 

As we did for the HBI, we study the saturation of the 
MTI using 2D and 3D simulations spanning a range of domain 
sizes L/H. Table |2] summarizes the simulations presented in 
this section. 

Since the linear dispersion relation successfully describes 
the nonlinear evolution and saturation of the HBI, it is a 
good place to begin our discussion of the MTI. The MTI is de- 
scribed by equation ([7| when dT/dz < 0. The linear evolution 
of the MTI is the opposite of that of the HBI: the MTI oper- 
ates when the temperature decreases with height, its fastest 
growing modes are the ones with wave vectors k parallel to 6, 
and the force that destabilizes the MTI is exactly that which 
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Figure 5. Evolution of the MTI with an initially horizontal magnetic field in a local, 2D simulation (simulation ml in Table [2I. Gray- 
horizontal lines show the transition to the buoyantly neutral layers described in § |3.2[ the color scale is identical to that in Figure [2] Initial 
perturbations grow by the mechanism described in § |2.2| (Fig.[l|>; rising and sinking plumes rake out the field lines until, by t = Stbuoy? they 
are mostly vertical. This configuration is, however, nonlinearly unstable to horizontal displacements, which generate a horizontal magnetic 
field and thus continually seed the MTI (see Fig.[6|. The result is vigorous, sustained convection in marked contrast to the saturation of the 
HBI in Fig. [2] In this local simulation, buoyant plumes accelerate until they reach the neutrally stable layers. The boundaries prematurely 
stop the growth of the MTI, and the local simulation under predicts the kinetic energy generated by the MTI (see Fig. [?]). 




Figure 6. Evolution of the MTI in a plasma with an initially vertical magnetic field in a local, 2D simulation (simulation m2 in Table |2j. 
This configuration is linearly stable according to equation ([7|, but there are zero-frequency, kx = 0, modes which do not have a restoring 
force. Physically, these correspond to horizontal motions which do not feel gravity/ buoyancy. Because of this zero frequency mode, small 
random initial perturbations add a horizontal component to the magnetic field, eventually rendering the plasma unstable to the MTI. This 
creates a feedback loop, allowing the MTI to generate vigorous, sustained convection; the HBI does not have this same feedback loop and 
so does not generate sustained turbulence (Fig. [2|. At late times, the results of this simulation with an initially vertical magnetic field are 
very similar to Figure [s] which starts with a horizontal magnetic field. 



stabilizes the HBI in its saturated state. Equation ([7|) shows 
that the maximum growth rate of the MTI goes to zero when 
bz = 1. By analogy with the HBI, it thus seems reasonable 
to expect that the MTI also saturates quiescently, by making 
the field lines vertical. 

The first three panels of Figure [5] show that this is 
nearly what happens. As the perturbations grow exponen- 
tially, the buoyantly rising and sinking blobs rake out the 
field lines, making them largely vertical. The growth rate of 
the MTI goes to zero when the field lines become vertical; 



since the velocities are still small at this point in the evolu- 
tion (^ 10~^Cs), one might expect the MTI to operate like 
the HBI and quiescently settle into this stable equilibrium 
state. Instead, however, the MTI drives sustained turbulence 
for as long as the temperature gradient persists. The plasma 
never becomes buoyantly stable, and the magnetic field and 
fluid velocities are nearly isotropic at late times. 

We can understand this evolution using the same ap- 
proach we employed for the HBI. Although the plasma in 
our MTI simulations never reaches a state in which the MTI 
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Table 2. Parameters for the MTI simulations (§|4.2| 



Name 


D 


res 


L/H 




Field Configuration 


ml 


2 


64 


0.033 


7.07 


horizontal 


m2 


2 


64 


0.033 


7.07 


vertical 


m3 


3 


64 


0.033 


7.07 


horizontal 


m4 


3 


64 


0.033 


7.07 


vertical 


m5* 


3 


128 


0.500 


0.31 


horizontal 


m6* 


3 


128 


1.400 


0.35 


horizontal 



The definitions of L, D, and n are the same as in Table ^ All 
simulations use the local setup (eq. [oj, except for the one with 
L/H = 1.4, which is global (eq. |10|. Each of these simulations was 
initialized with a weak magnetic field B /^/Ait = 10~^ with the ori- 
entation indicated in the table. *We also repeated simulations m5 
and m6 with initial field strengths B j \fA^ = 10"^, 0.0014, 0.0245 



growth rate is zero, examining the properties of this state is 
very instructive. The equilibrium state of the MTI with hz — 1 
(i.e., a vertical field) has precisely the same dispersion rela- 
tion as the saturated state of the HBI, given by equation ( 12 ). 



There are again zero frequency (neutrally stable) modes of the 
dispersion relation which correspond to horizontal perturba- 
tions to the equilibrium state of the MTI; these experience no 
restoring force, because the restoring force is buoyant in na- 
ture and unaffected by horizontal displacements [j Critically, 
however, these zero frequency perturbations now add a hor- 
izontal component to the magnetic field, pulling the plasma 
out of the equilibrium state and rendering it unstable to the 
MTI. 

Figure [G] vividly illustrates this process. This figure shows 
a simulation that starts with the linearly stable equilibrium 
state of the MTI (a vertical magnetic field), seeded with the 
same highly subsonic velocity perturbations as before. The 
compressive component of the perturbation rapidly damps 
because the Mach number is small, and buoyancy traps the 
vertical components of the perturbation in small-amplitude 
oscillations, as predicted by equation (12). The incompres- 



sive, horizontal displacements propagate freely, however, and 
by t = 7.5tbuoy, they have noticeably changed the local orien- 
tation of the magnetic field. The plasma is no longer in its sta- 
ble equilibrium state, and by t = 20tbuoy, it is clear that this 
process has excited the MTI. At late times, the simulations 
initialized with horizontal (linearly unstable; Fig.|5]) and ver- 
tical (linearly stable; Fig. [6| magnetic fields are qualitatively 
indistinguishable. Thus, although equation ([7| shows that a 



^ This conclusion is more subtle than the analogous argument for 
the HBI, because the equilibrium state has a nonzero heat fiux. 
If the horizontal displacements aren't incompressive, the field lines 
could pinch together, heating and destabilizing parts of the plasma. 
These perturbations do not appear in eq \12\ because we have 
taken the Boussinesq limit. It is in principle possible that such 
compressive perturbations contribute to destabilizing the bz = 1 
MTI state in our simulations, but all of our analysis is consistent 
with the neutrally stable zero frequency perturbations being the 
critical ingredient. 



plasma with dT/dz < is linearly stable if the magnetic field 
is vertical, that configuration is nonlinearly unstable. 

The nonlinear instability of the bz — 1 state of the 
MTI precludes the magnetic saturation channel. The zero 
frequency horizontal motions generate a horizontal magnetic 
field component from the vertical magnetic field, seeding the 
instability and closing the dynamo loop. This continuously 
drives the MTI and generates sustained turbulence. Without 
a linear means to saturate, the MTI grows until nonlinear ef- 
fects can compete with the linear instability, which requires 

V ^ Cs. 

Figure [T] shows the Mach number in 3D MTI simulations 
as a function of time, for different sizes of the computational 
domain L/H. The MTI buoyantly accelerates rising or sink- 
ing fluid elements. For simulations with L <^ the velocities 
generated by the MTI are artiflcially suppressed because the 
small size of the computational domain prematurely stops the 
buoyant acceleration (Fig. [5]); the results in this Figure [t] are 
reasonably consistent with mixing length estimate v ^ \fTj. 
By contrast, for simulations with L > H, the Mach numbers 
approach ^ 1 so that the MTI taps into the full buoyant force 
associated with the unstable temperature gradient. This is 
only true, of course, because our boundary conditions fix the 
temperature at the top and bottom of the computational do- 
main. If the temperatures were free to vary, the MTI could 
saturate by making the plasma isothermal. Which of these 
saturation mechanisms is realized in a given astrophysical sys- 
tem will depend on the heating and cooling mechanisms that 
regulate the temperature profile of the plasma. 

Figure [S] shows the evolution of the magnetic field ori- 
entation (left panel) and energy (right panel) as a function 
of time in our 3D MTI simulations (for two different val- 
ues of L/H). The sustained, vigorous turbulence generated 
by the MTI in the nonlinear regime (t > lOtbuoy) rapidly 
amplifies the magnetic field. The magnetic and kinetic ener- 
gies reach approximate equipartition in our simulations, with 
B^/Stt ^ 0.1 pv^. The kinetic and magnetic energies gener- 
ated by the MTI together contribute ^ 5-10% of the pressure 
support in its saturated state. 

The evolution of the magnetic field geometry, shown in 
the left panel of Figure [8] nicely illustrates the transition of 
the MTI from the linear to nonlinear regime. During the linear 
phase of the instability (t < lOtbuoy), the plasma accelerates 
toward the equilibrium state with 6^ = 1- After the evolution 
becomes nonlinear (t > lOtbuoy), however, the MTI drives 
sustained turbulence, which nearly isotropizes the magnetic 
field. 

While the HBI works to insulate the plasma against ver- 
tical energy transport, the MTI enhances it. Figure [S] shows 
that the conductive flux through the plasma is slightly greater 
than ^ 1/3 of the fleld free value KqVT (because (bl) ^ 0.4). 
Moreover, the MTI leads to large fluid velocities and corre- 
lated temperature and velocity perturbations — hot pockets 
of plasma rise, while cool pockets sink. These imply that the 
MTI drives an efficient outwards convective heat flux. Fig- 
ure [o] shows that this flux can be ^ 1.5% of pcf, consistent 
with the Mach numbers of 0.2 in Figure [7| This convective 
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Figure 7. Volume-averaged Mach numbers of the turbulence gen- 
erated by the MTI in 3D simulations, for three different values of 
the size of simulation domain relative to the temperature scale- 
height L/H. The MTI buoyantly accelerates the unstable fluid 
elements; if the size of the simulation domain is smaller than a 
scale-height, the boundaries suppress the growth of the instabil- 
ity. Local simulations with L/H = 1/30 therefore strongly under 
predict the strength of the turbulence generated by the MTI. In 
global simulations with L > H, the MTI leads to turbulence with 
average Mach numbers of ~ 0.2; the velocity distribution extends 
up to ~ 5 times the mean. 



flux is probably not large enough to influence the thermody- 
namics of the ICM, but it could be important in other envi- 
ronments where the MTI can operate, such as the interiors of 
white dwarfs and neutron stars ( Chang et al.||2010[ ). 



5 INTERACTION WITH OTHER SOURCES OF 
TURBULENCE 

Thus far, we have described the development and saturation 
of buoyancy instabilities only in the idealized case of an oth- 
erwise quiescent plasma. In a more astrophysically realistic 
scenario, however, other processes and sources of turbulence 
may also act on the plasma and the resulting dynamics can 
be more complicated. For example, the evolution of the HBI 
won't simply proceed until the growth rate is everywhere zero. 
Instead, we expect the saturated state to involve a statisti- 
cal balance among the various forces; this balance depends 
on the buoyant properties of the plasma and provides a test 
of our understanding of the nonlinear behavior of the HBI 
and MTI. Furthermore, any change in the saturated state of 
the plasma due to the interaction between the HBI/MTI and 
other sources of turbulence could change the astrophysical 
implications of these instabilities. 
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Figure 8. Left panel: Evolution of the magnetic field orientation 
in 3D MTI simulations, for two different values of the size of simula- 
tion domain relative to the temperature scale-height L/H. During 
the linear phase of evolution {t < lO^buoy)? the MTI drives the 
plasma toward a nearly vertical magnetic field, i.e., 6^ ~ 1. When 
the instability becomes nonlinear, however, (near t ~ lOtbuoy) the 
evolution changes. Unlike the HBI, the plasma never settles into 
an equilibrium state; instead the MTI drives vigorous turbulence. 
This turbulence amplifies and nearly isotropizes the magnetic field. 
Right panel: Evolution of the magnetic (dashed line) and kinetic 
(solid line) energy in local (L/H = 1/2) 3D MTI simulations, with 
different initial field strengths. These local simulations have a pos- 
itive entropy gradient, but under-predict the magnetic and kinetic 
energy produced by the MTI. The magnetic energy in the saturated 
state approaches ~ 10% x p-u^. 



We choose to explore the interaction between buoyancy 
instabilities and other sources of turbulence using the ide- 
alized, isotropic turbulence model described in section |3.4[ 
While this model glosses over the details of what generates 
the turbulence, we hope that it captures the essential physics 
of the problem, allowing us to study the effect of turbulence 
without unnecessarily restricting our analysis to specific ap- 
plications. We intend to specialize to specific sources of turbu- 
lence in future work, but our present analysis should apply in 
the ICM, accretion disks, and anywhere else the assumptions 
summarized in section [2?2] apply. 

In order to characterize the turbulence, we define a 
timescale for it to influence the plasma. We define this "distor- 
tion time" in terms of the spatial velocity spectrum: tdist(^) = 
i/Sv{i), where 



6v(i) 



J {Sv{k)Y S{\k\-27T/i) 



(fk 



(27r)3 



1/2 



(13) 



and 5v{k) is the Fourier transform of the velocity field. We 
expect the relevant parameter describing the importance of 
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Figure 9. Convective energy fluxes generated by the MTI. The 
large turbulent velocities associated with the MTI (Fig. [7| lead to 
eflicient convective transport of energy, at a reasonable fraction of 
the maximal value, pc^. In these MTI simulations, the conductive 
energy flux is ~ 0.4 times the fleld-free value and always dominates 
the convective energy flux for galaxy cluster conditions. 



the turbulence to be the ratio of the timescales tbuoy/tdist- 
This represents a dimensionless strength of the turbulence; if 
^buoy/tdist ^ 1, the turbulence displaces fluid elements faster 
than buoyancy can restore them, and we expect the velocities 
and magnetic field lines to become isotropic. In the opposite 
limit, buoyancy still plays an important role in the evolution 
of the plasma. The scale dependence of the distortion time 
makes the ratio tbuoy/^dist a function of scale. We define this 
ratio at the scale where the velocity spectrum of the injected 
turbulence peaks. This is roughly consistent with the driving 
scale of the turbulence and typically represents the scale with 
the most energy. 

In the following sections, we study the transition from a 
state dominated by buoyancy to one dominated by isotropic 
turbulence using a number of simulations of the HBI and 
MTI, with turbulence in the range 0.1 < tbuoy/tdist ^10. 



5.1 Effect of Turbulence on the HBI 

Figure [To] shows representative snapshots of the temperature 
and magnetic field lines in saturated states of our HBI simula- 
tions with turbulence; the strength of the injected turbulence 
increases from the top left panel through the bottom right 
(the labels correspond to the values of tbuoy/^dist)- When tdist 
is long compared to the buoyancy time, as in the first panel 
of Figure [To] the turbulence is weak; the HBI therefore dom- 
inates and the evolution of the plasma is similar to that de- 
scribed in section 14.11 This saturated state of the HBI feels 



Table 3. Parameter study for the HBI simulations with turbulence 



D (2) res (64) L (0.1) H (2.0) Bq (10"^) ko 



1.0 
3.0 



128 



256 



1.0 
0.3 



1.0 
1.0 



10-3 
3 X 10-4 
3 X 10-4 



The simulations were initialized with our local setup (eq.[8| and an 
initial magnetic fleld strength Bq. Each simulation was performed 
on uniform Cartesian grids of side L, resolution res and dimension 
D. We varied the size of the simulation domain L (scaling the con- 
ductivity as described in § |3.3[ ), the plasma scale height H and the 
initial magnetic fleld strength Bq; the flducial values for these pa- 
rameters are included in the table header ( — indicates the flducial 
value). For each entry in the table, we performed simulations with 
both initially horizontal and vertical magnetic flelds. As described 
in the § |5.1[ these simulations include isotropic turbulence injected 
at the scale ko = 27t/L x ko and with a range of turbulent energy 
injection rates to give 0.1 < ^buoy/^dist ^ 10- 



a buoyant restoring force which resists vertical displacements 
with a force per unit mass /buoy = ^buoy-f^- As we increase 
the strength of the applied turbulence in the following panels 
of Figure |10| the vertical displacements grow, and the field 
deviates more strongly from the bz = equilibrium state of 
the HBI. When tdist is short compared to the buoyancy time, 
as in the last panel of Figure |10| turbulence can displace the 
fluid elements faster than buoyancy can restore them. The 
turbulence then dominates the evolution of the plasma, tan- 
gling and isotropizing the field lines. 

Figure [To] also shows the length-scale dependence of the 
transition from an HBI to a turbulence dominated state. It is 
clear in the second and third panels that the HBI has globally 
rearranged the field lines, but that the turbulence is increas- 
ingly efficient at smaller scales. This is a consequence of the 
fact that turbulence typically perturbs the plasma in a scale- 
dependent way, while the buoyant restoring force of the HBI 
does not. If the turbulence follows a Kolmogorov cascade, the 
force uj'^/k oc k^^^ increases with decreasing scale, so the 
turbulence will always win on sufficiently small length-scales. 
It is therefore somewhat ambiguous whether turbulence or 
the HBI dominates a certain configuration, as the answer will 
typically depend on scale. As mentioned earlier, we skirt this 
issue by defining tdist at the scale where the velocity spectrum 
of the injected turbulence peaks. When assessing the astro- 
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Figure 10. Snapshots of the saturated states of our 2D HBI simu- 
lations with externahy driven turbulence. Colors show the temper- 
ature (increasing from blue to red), and black lines show magnetic 
field lines. Each panel is labeled with the dimensionless strength of 
the turbulence in the simulation, tbuoy/^dist? defined in §[5] When 
^buoy ^ ^dist5 cis in the top two panels, the HBI dominates the 
evolution of the plasma. When tbuoy ^ ^dist? the turbulence can 
isotropize the magnetic field, but it does so in a scale-dependent 
way with the large scales retaining memory of the horizontal field 
imposed by the HBI. 



physical importance of the HBI, it is important to keep this 
scale in mind. If the scale where the turbulent energy spec- 
trum peaks is smaller than the temperature gradient length 
scale, the HBI may still insulate the plasma against conduc- 
tion, even if the field lines are isotropized on smaller scales. 

In order to quantify the transition from an HBI- 
dominated configuration to one dominated by turbulence, we 
measure the mean orientation of the magnetic field via the 
volume average of b^. The saturated value for this quantity 
approaches zero when the HBI dominates, and 1/D when the 
magnetic field is isotropic, where D is the number of dimen- 
sions in the simulation. 

Figure [TT] shows the saturated field angle as a function 
^buoy/tdist- The points in this figure represent simulations with 
different driving scales ko , different dimensionality, and differ- 
ent buoyancy times tbuoy (Table [s] summarizes our parameter 
study). Symmetry of the coordinate axes requires that bl = 
1/2 in 2D or 1/3 in 3D if the magnetic field is isotropic. To 
include both our 2D and 3D simulations on the same plot, we 
shift our 3D values of (hi) by a factor of 3/2. To within the 
scatter shown in Figure [TT] we find that the saturated value of 
(bl) depends only on the ratio tbuoy/tdist: for tbuoy ^ tdist the 
turbulence is strong and the field becomes relatively isotropic 



while for tbu 



tdist the isotropic turbulence is weak and 
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Figure 11. Saturated magnetic field orientation as a function of 
the strength of the externally driven turbulence, tbuoy/^dist? for 



local HBI unstable atmospheres with tbu 



1, and v^. The 



the HBI drives the magnetic field to become relatively hori- 



thick gray line is an isotropic magnetic field in 2D. Colored points 
represent simulations with turbulence driven on different scales. 
Squares mark 3D calculations; the values of (6^) for the 3D sim- 
ulations have been shifted by a factor 3/2 since isotropy implies 
6^ = 1/2 in 2D but 6^ = 1/3 in 3D. Error bars represent la statis- 
tical fluctuations in bz and tdist- 



zontal. The fact that the transition between these two states 
occurs around tbuoy ^ tdist suggests that our definition of 
tdist, though somewhat arbitrary, is reasonable. 

The bulk of the simulations in Figure [TT] are 2D, and 
we do not have any 3D simulations in the very weak turbu- 
lence limit. These simulations are computationally expensive, 
both because of conduction and because we have to run for 
a long time for turbulence and the HBI to reach a statisti- 
cal steady state; using 2D simulations allowed us to explore 
a larger fraction of the interesting parameter space. While 
the development of turbulence is very different in two and 
three dimensions, the HBI is essentially two-dimensional in 
nature. Moreover, the key dynamics governing the interaction 
between the HBI and the turbulence are dominated by the 
energy-containing scale of the turbulence — the precise power- 
spectrum of the fluctuations (which differs in 2D and 3D) is 
less critical. Scaling for dimension, we find that the saturated 
states of our 2D and 3D simulations are nearly identical. We 
thus believe that results in Figure [TT] in the weak turbulence 
limit are a good description of the magnetic field structure in 
3D systems as well. 

These results on the interaction between the HBI and 
other sources of turbulence support an analogy between the 
saturated state of the HBI and ordinary, adiabatic stable 
stratification. The most important effect of the HBI in the 
saturated states of our simulations is to inhibit mixing in the 
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Figure 12. Magnetic field orientation as a function of time in 2D 
simulations with turbulence and the HBI for different initial mag- 
netic field orientations; tbuoy/^dist = 0-8- Left panel: Simulations 
with anisotropic thermal conduction. Right panel: Adiabatic simu- 
lations with no conduction (and thus no HBI). Different curves rep- 
resent simulations with initially vertical and horizontal magnetic 
field lines, respectively. Conducting simulations eventually reach 
the same saturated state independent of the initial magnetic field 
direction. Adiabatic simulations are very similar to the conducting 
ones if the field is initially horizontal, highlighting the fact that in 
the saturated state of the HBI the plasma is buoyantly stable and 
behaves dynamically like an adiabatic fiuid. 



direction of gravity. Ordinary stable stratification also inhibits 



vertical mixing and, as suggested by | Sharma et al. (2009b), 
our parameter tbuoy/tdist is analogous to the Richardson num- 
ber used in the hydrodynamics literature. To further expand 
on this analogy. Figure [12] shows a comparison of anisotropi- 
cally conducting (left panel) and adiabatic (right panel) sim- 
ulations with equal values of tbuoy/tdist- For the adiabatic 
simulations, we define tbuoy = ^ad, i-e., using the entropy gra- 
dient rather than the temperature gradient. The left panel 
of Figure [12] shows simulations with the same injected turbu- 
lence, but initialized with vertical or horizontal magnetic field 
lines. In these simulations, the saturated state is independent 
of the initial condition. That is, the interaction between tur- 
bulence and the HBI leads to a well defined magnetic field 
orientation that is independent of the initial field direction. 

By contrast, in the adiabatic simulations (right panel of 
Fig. 12), the final magnetic field orientation depends on the 
initial field direction. For an initially vertical field in an adia- 
batic plasma, the turbulence slowly isotropizes the magnetic 
field direction. However, the adiabatic simulations with ini- 
tially horizontal field lines reach a saturated state that is very 
similar to that of the HBI simulations. In the adiabatic simu- 



lations, the magnetic field is essentially passive, but it traces 
the fluid displacements. The stable stratification competes 
with the turbulence and sets a typical scale for vertical dis- 
placements in the saturated state. This scale, in turn, deter- 
mines the magnetic field geometry. The magnetic field plays 
no dynamical role in this process. The fact that the anisotrop- 
ically conducting simulations reach the same statistical steady 
state highlights that the saturated state of the HBI behaves 
dynamically very much like an adiabatic, stably stratified, 
plasma. 



5.2 Effect of Turbulence on the MTI 

To complete our analysis, we study how externally imposed 
isotropic turbulence affects the saturation of the MTI. Fig- 
ure [13] shows the volume averaged magnetic field orientation 
as a function of time in MTI simulations with additional tur- 
bulence, for different values of the strength of the turbulence 
tbuoy /tdist- These are local simulations (with domain sizes 
L/H — 0.5; initialized using eq. (|9|) which have the peda- 
gogical advantage of a positive entropy gradient, but which 
under-predict the kinetic energy generated by the MTI. We 
drive the turbulence at relatively small scales {kL/2Ti — 8) 
so that subsonic turbulence can still satisfy tbuoy /tdist > 1- 
Both our driven turbulence and the MTI tend to isotropize 
the magnetic field, so it is not a priori clear whether the field 
orientation is a good indication of the importance of turbu- 
lence relative to the MTI. Although Figure [13] shows that 
there is no strong dependence of the saturated field orienta- 
tion on the strength of the turbulence, the time dependence of 
the field orientation clearly shows the effects of the turbulence 
on the MTI. 

In general, the evolution of the MTI proceeds through 
two stages: there is a linear phase, where the plasma acceler- 
ates toward its nominal stable state (which has a vertical mag- 
netic field), and a nonlinear transition to the saturated state, 
where the strong turbulence generated by the MTI isotropizes 
the velocities and field lines. In the absence of additional tur- 
bulence, the linear phase is characterized by field lines that 
are primarily in the direction of gravity (Fig. [s]). Figure 13 
shows that additional sources of strong (rapidly shearing) tur- 
bulence suppress this linear phase of the MTI. Indeed, the 
evolution of the field angle with time in our strongest turbu- 
lence simulations (tbuoy /tdist = 4.7) is quite similar to what 
we find in simulations of an adiabatic plasma in which the 
MTI is not present. 

It would, however, be incorrect to conclude from Fig- 
ure [13] that the MTI is unimportant if there are other strong 
sources of turbulence in the plasma. The fundamental rea- 
son for this is that the growth of the MTI does not depend 
significantly on scale, while the effects of the other sources 
of turbulence do. Figure [M] shows velocity spectra for the 
simulations in Figure [Ts] for comparison we also show the ve- 
locity spectra in adiabatic simulations, which correspond to 
the power spectra produced solely by the injected turbulence. 
Figure ^] demonstrates that even in simulations with very 
strong imposed turbulence (tbuoy /tdist = 4.7) there is still 
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Figure 13. Magnetic field orientation as a function of time in 3D 
simulations with turbulence and the MTI, for different values of the 
strength of the turbulence tbuoy/^dist- In all of the simulations, the 
magnetic field is relatively isotropic in the saturated state. How- 
ever, the early-time 'overshoot' towards a vertical magnetic field 
due to the MTI is suppressed in the presence of strong turbulence 
with tbuoy ^ ^dist- Thcsc relatively local simulations (L = 0.5if) 
underestimate the kinetic energy generated by the MTI (Fig. [?]) 
and therefore the strength of the turbulence required to infiuence 
it. 



Figure 14. Velocity spectra in simulations with turbulence and the 
MTI (solid lines), for different values of the strength of the injected 
turbulence tbuoy/^dist- The turbulence is driven at kL /2tt = 8. For 
comparison, dashed lines show the power spectra in adiabatic simu- 
lations, in which the injected turbulence is the only source of power. 
Even when tbuoy/^dist ^ 1 at the driving scale of the turbulence, 
so that one might expect the MTI to be suppressed, the MTI pro- 
duces significant turbulent energy on larger scales kL /2t: ~ 2 — 10. 
Suppressing the MTI on all scales in a plasma would thus require 
^buoy/^dist ^ 1 on all scalcs smaller than the scale-height. 



significant excess power on the largest scales in the compu- 
tational domain {kL /2ti < 10). This large-scale power is due 
to the MTI. Moreover, the turbulent energy due to the MTI 
dominates the total turbulent kinetic energy in the plasma. 
These results highlight that 'strong' turbulence is a scale- 
dependent statement. Suppressing the MTI requires having 
^buoy/^dist ^ 1 ou all scalcs, up to the temperature/pressure 
scale-height of the plasma. Because the MTI itself generates 
nearly sonic velocities, this suppression would require close 
to supersonic turbulence. In practice, it is therefore unlikely 
that additional sources of turbulence can fully suppress the 
MTI in most astrophysical environments where it is likely to 
occur (e.g., accretion disks and galaxy clusters). 



6 DISCUSSION 

The motion of electrons along, but not across, magnetic 
field lines in dilute, magnetized plasmas produces efficient, 
anisotropic transport of heat. Such plasmas are therefore non- 
adiabatic, and the standard analysis of buoyancy (or con- 
vective) instabilities does not necessarily apply. Quantita- 
tively, conduction plays an essential role on scales less than 

7 (Aif)^/^, where A is the electron mean free path and H is the 
plasma scale height. In this "rapid conduction limit," the tem- 
perature gradient, rather than the entropy gradient, dictates 



the stability of the plasma, and the plasma is unstable for ei- 
ther sign of the temperature gradient ( |Balbus|2000 Quataert 
2008| ) . The convective instability in this limit is known as the 



HBI (MTI) when the temperature increases (decreases) with 
height. 



Parrish k Stone (2005 2007) and Parrish k Quataert 



( 2008 ) extended the original linear analysis of the MTI and 
HBI into the nonlinear regime using numerical simulations. 
In this paper, we have reconsidered the nonlinear saturation 
of the HBI and MTI. Our work adds to previous investiga- 
tions because we have identified a key difference between the 
two instabilities and are able to understand the nonlinear 
behavior of the MTI more completely. This paper therefore 
represents a significant change in our understanding of the 
possible astrophysical implications of the MTI (but not the 
HBI). We have also studied the effect of an external source 
of turbulence on both the MTI and HBI. We conclude that 
other sources of turbulence in a plasma can change the satu- 
ration of the HBI, but that it is much harder to disrupt the 
MTI. Below we summarize our results and discuss their as- 
trophysical implications, focusing on the intracluster medium 
(ICM) of galaxy clusters. 
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6.1 HBI 

The HBI occurs whenever the temperature increases with 
height in an anisotropically conducting plasma. Plasmas in 
the rapid conduction limit are in general linearly unstable, 
unless the magnetic field lines are horizontal (eq.[7|. Horizon- 
tal field lines represent a fixed point in the evolution of the 
plasma: if the plasma somehow reaches such a state (and is 
not perturbed away by another process), it will remain there 
forever. A horizontal magnetic field is therefore a natural sat- 
uration channel for the HBI. 

This "magnetic" saturation mechanism can be under- 



than its age; the HBI removes thermal conduction as a source 
of energy for the cores, potentially exacerbating the cooling 



stood using the linear dispersion relation of the plasma (§ 4.1 ). 
Starting from a linearly unstable state, the HBI induces both 
horizontal and vertical motions in the plasma. As the HBI 
develops, however, the vertical motions become trapped in 
internal gravity waves. These waves decay, leaving only the 
horizontal motions at late times; thus, the fluid velocities are 
very anisotropic in the saturated state. Since the horizontal 
motions don't incur a buoyant response, the horizontal dis- 
placements can be very large. These motions stretch out the 
magnetic field lines, amplifying and reorienting them, and 
drive the plasma towards its stable equilibrium with horizon- 
tal magnetic field lines. These horizontal motions therefore 
drive the nonlinear evolution and saturation of the HBI. 

Since the saturation of the HBI is dominated by hori- 
zontal displacements, the key dynamics all occur at approxi- 
mately the same height in the atmosphere. The saturation of 
the HBI is thus essentially local in nature. We have demon- 
strated this explicitly by carrying out simulations with dif- 
ferent domain sizes relative to the plasma scale- height; the 
results of these simulations are very similar (Fig. [4|. 

The growth rate of the HBI decreases dramatically as 
the instability progresses. The HBI therefore saturates qui- 
escently, and the velocities in the saturated state are very 
subsonic (Fig.[3|. The saturation is driven by horizontal mo- 
tions with nearly constant velocities, so the nonlinear mag- 
netic field amplification is approximately linear, rather than 
exponential, with time. These findings are consistent with 
those of'Parrish & Quataert' (2008). We note that the HBI is 
very much unlike adiabatic convection, which can only satu- 
rate by changing the thermal state of the plasma and therefore 
generates vigorous turbulence. The key difference between the 
HBI and adiabatic convection is that the source of free energy 
for the HBI is a conductive heat flux through the plasma, not 
merely the existence of a temperature gradient. Since the heat 
flux can be suppressed by rearranging the magnetic field, the 
HBI has a magnetic saturation channel that is not available 
to adiabatic convection. 

The astrophysical implications of the HBI follow immedi- 
ately from the nature of its saturated state. By reorienting the 
magnetic field lines, the HBI dramatically reduces the conduc- 
tive heat flux through the plasma. The HBI should operate 
in the innermost ^100-200 kpc in the intracluster medium 
of cool-core galaxy clusters, where the observed temperature 
increases outward. As noted in iParrish & Quataert (2008), 



flow problem ( Parrish et al.|2009 ). 

Our results demonstrate that the saturated state of the 
HBI is buoyantly stable. This may not seem surprising, be- 
cause it is exactly what one would expect if the ICM were 
adiabatic. The ICM is not, however, adiabatic, and thermal 
conduction would render it buoyantly neutral to vertical dis- 
placements if the magnetic field lines were tangled. The satu- 
rated state of the HBI is buoyantly stable only because of the 
nearly horizontal magnetic field lines (that are perpendicular 
to the temperature gradient). The HBI therefore inhibits ver- 
tical mixing and allows for the existence of weakly damped 
internal gravity waves in the plasma. 



Sharma et al. (2009b) noted that the stable stratifica- 



tion associated with the saturated state of the HBI competes 
with other sources of turbulence in a well-defined way. This 
competition can be understood using a modified Richardson 
number tbuoy/tdist, where tbuoy is the timescale for the HBI to 
grow and tdist is a characteristic "distortion time," or "eddy 
turnover time," of the turbulence. When tbuoy ^ tdist, the 
turbulence can isotropize the plasma and remove all traces 
of the HBI (Fig. 11). When tbuoy ^ tdist, the saturated state 
of the plasma represents a statistical balance between turbu- 
lence and the HBI, with the magnetic field becoming more 
horizontal, and the plasma more HBI-dominated, for smaller 
values of tbuoy /tdist- The strength of other sources of turbu- 
lence is therefore crucial for understanding the astrophysical 
implications of the HBI. 

Figure [TT] provides a very simple mapping between the 
properties of the turbulence and the magnetic field geome- 
try in the plasma. Given the strength of the turbulence and 
the thermal state of the plasma, this figure provides a recipe 
for determining the mean geometry of the magnetic field, and 
therefore the effective conductivity of the plasma. This can be 
used to interpret observational results or to construct semi- 
analytic models of anisotropic conduction for use in cosmo- 
logical simulations. 

Turbulence in the ICM is currently poorly constrained, 
and thus it is difficult to determine precisely how important 
the HBI is for the evolution of clusters. Reasonable estimates 
suggest that tbuoy /tdist ^ 1, but more detailed simulations 
of clusters are required to determine this ratio more pre- 
cisely. Future observations of the ICM with space-based x-ray 
calorimeters will place observational constraints on the level 
of turbulence. In addition, Faraday rotation measurements of 
the ICM will measure the orientation of the magnetic fields in 
clusters and constrain the role of the HBI (Bogdanovic et al.| 
2010) (Pfrommer & D ursi (2010) describe another mechanism 
to measure the magnetic orientation in the ICM). Even if the 
turbulence is strong (tbuoy /tdist ^ 1), the driving scale and 
filling factor of the turbulence may allow the HBI to dom- 



inate on some scales or at some locations (see Fig. 10 and 
associated discussion in §|5.1|. As suggested by Parrish et al. 



(2010) and Ruszkowski & Ohl ( 2010 ) , the interaction between 



this is precisely where the cooling time of the ICM is shorter 



turbulence and the HBI might be part of a feedback loop for 
the thermal evolution of the ICM. 
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6.2 MTI 

The nonlinear evolution of the MTI is more complex than 
that of the HBI. Just hke the HBI, the MTI has linearly stable 
equilibria, but they are transposed: the linearly stable equi- 
librium states of the MTI have vertical field lines. We have 
shown, however, that the linearly stable equilibrium states 
turn out to be nonlinearly unstable; i.e., they are unstable to 
perturbations with a finite amplitude. This nonlinear instabil- 
ity arises because neutrally buoyant, horizontal displacements 
add a horizontal component to the magnetic field. This takes 
the plasma out of its linearly stable state and re-seeds the 
instability (Fig. [g] and § 4.2). As a result, the linearly stable 
equilibrium states of the MTI do not represent fixed points in 
the evolution of the instability, and the MTI cannot saturate 
simply by reorienting the magnetic field. 

This difference eliminates the quiescent, magnetic satura- 
tion channel for the MTI and dramatically changes its evolu- 
tion. Without a linear means to saturate, the instability grows 
until nonlinear effects dominate, which occurs when v ^ Cg. 
Unlike the HBI, the MTI therefore drives strong turbulence 
and operates as an efficient magnetic dynamo, much more 
akin to adiabatic convection. The astrophysical implications 
of the MTI are therefore entirely different from those of the 
HBI. The kinetic and magnetic energy generated by the MTI 
can contribute a significant (up to ten percent) non-thermal 
pressure support to the plasma in the saturated state. This 
is consistent with observational constraints on non-thermal 
pressure support in the ICM near the virial radius ( |George| 
et al.|2009| . This non-thermal pressure support may have con- 
sequences for mass estimates of clusters, which often rely on 
the assumption of hydrostatic equilibrium with thermal pres- 
sure support. Note, in particular, that the MTI is predicted 
to be present at precisely the same radii (> the scale radius) 
to which x-ray and SZ mass measurements are most sensitive. 

Because the MTI operates by buoyantly accelerating fluid 
elements until they approach the sound speed, the results of 
numerical simulations of the MTI are sensitive to the size of 
the computational domain. The boundaries of the domain can 
artificially suppress this acceleration, and simulations with 
sizes smaller than a scale height under-predict the kinetic en- 
ergy generated by the MTI (this was the case in the original 
MTI simulations of iParrish & Stone||2005| 120071). The non- 



linear development of the MTI is therefore quite sensitive to 
the global thermal state of the plasma, and an understanding 
of the MTI requires more careful numerical simulations than 
are needed for the HBI. 

The saturated state of the MTI corresponds to a largely 
isotropic magnetic field, with a slight but persistent vertical 
(or radial) bias; this bias is robust even in the presence of 



other sources of strong turbulence (Fig. 13). We find that the 



magnetic energy generated by the MTI saturates at about 
30% of the kinetic energy (Fig.[8|. However, it may be difficult 
to observationally distinguish the turbulence generated by the 
MTI from that generated by other processes. 

The large velocities generated by the MTI, along with 
correlations between the temperature and velocity perturba- 



tions, imply that the MTI drives a large convective heat flux, 
^ 1.5% X pel (Fig.§. While this convective flux is probably 
too small to alter the thermal evolution of the ICM, it could 
be important in higher density astrophysical plasmas, where 
the electron mean free path is smaller and conduction isn't 
as efficient. 

While the MTI cannot saturate by reorienting the mag- 
netic field, it can saturate by making the plasma isothermal. 
Simulations with Neumann boundary conditions in which the 
temperature at the boundaries is free to adjust find that this 
is the case; the atmosphere becomes isothermal before the 
MTI has a chance to develop ( [Parrish et aL]|2008t . We fixed 
the temperature at the top and bottom boundaries of our 
simulations; this is partially motivated by the fact that many 
galaxy clusters in the local universe are observed to have non- 
negligible temperature gradients. 



Sharma et al. ( 2008 ) carried out numerical simulations of 



the MTI in spherical accretion flows and found nearly radial 
magnetic fields, with modest turbulence. This quasi-linear 
saturation of the MTI might seem to contradict the results 
presented in this paper. Note, however, that in the Bondi in- 



flow studied by Sharma et al. (2008|, the plasma undergoes at 
most 5 — 10 MTI growth times before flowing in. After 10 
growth times, our simulations also show approximately radial 
field lines and modest turbulence (Fig.jsJ. Moreover, the sim- 
ulations of Sharma et al. ( 2008 ) covered a very large dynamic 



range in radius and may have lacked the resolution to see the 
full nonlinear development of the MTI. 

The MTI growth time in the outer parts of the ICM 
is about 1 Gyr. Although our typical MTI simulations take 
^10-20 growth times to saturate, this does not preclude the 
importance of the MTI in galaxy clusters. Figure [S] shows 
that there is a long, linear ramp-up phase where the instabil- 
ity grows from the tiny perturbations we apply to the non- 
linear state. Astrophysical perturbations are unlikely to be 
this subsonic. Figure [13] shows that the MTI can saturate in 
^2-5 growth times when subjected to larger perturbations, 
suggesting that the MTI is likely to become nonlinear in the 
outer parts of clusters. Cosmological simulations will be re- 
quired to fully understand the implications of the MTI for 
galaxy clusters. 
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